//------------------------------------------------------------------
//  POLY.c
//------------------------------------------------------------------
#include <stdio.h>
//-------------------------------------------------------------------
#include "Top.h"
#include "Bigint.h"
#include "BigintEC.h"
#include "Poly.h"
#include "PolyEC.h"
//------------------------------------------------------------------
//  Define irreducible polynomial for field here.  Use table, or evaluate your
//	own, but the math won't work if it' isn't really a prime polynomial.  
//	NOTE: the value you pick must be consistent with NUMBITS.
//------------------------------------------------------------------
FIELD2N poly_prime = {0x00000008,0x00000000,0x00000000,0x00000000,0x00000000,0x000000c9}; /*163*/
//FIELD2N poly_prime = {0x08000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x00000000,0x000010a1}; /*283*/
/*FIELD2N poly_prime = {0x08000000,0x00000000,0x00000000,0x00000000,0x000000B1};*/  /*155*/
/*FIELD2N poly_prime = {0x08000000,0x00000000,0x00000000,0x40000000,0x00000001};*/  /*155*/
/*FIELD2N poly_prime = {0x00000040,0x00000000,0x00000000,0x02000000,0x00000001};*/	/*134*/
/*FIELD2N poly_prime = {0x00000008,0x00000000,0x00000000,0x00000000,0x0000010D};*/	/*131*/
/*FIELD2N poly_prime = {0x00000004,0x00000000,0x00000000,0x00000000,0x00000009};*/	/*130*/
/*FIELD2N poly_prime = {0x00000002,0x00000000,0x00000000,0x00000000,0x00000021};*/	/*129*/
/*FIELD2N poly_prime = {0x00000001,0x00000000,0x00000000,0x00000000,0x00000087};*/	/*128*/
//FIELD2N poly_prime = {0x00008000,0x00000000,0x00000000,0x00000401};*/				/*111*/
/*FIELD2N poly_prime = {0x20000000,0x00000000,0x00000005};*/	/*93*/
/*FIELD2N poly_prime = {0x0000800,0x00000000,0x0000004b};*/		/*75*/
/*FIELD2N poly_prime = {0x0000020,0x00000000,0x00000065};*/		/*69*/
/*FIELD2N poly_prime = {0x00000002,0x00000000,0x00000001b};*/	/*65*/
/*FIELD2N poly_prime = {0x00000001,0x00000000,0x00000001b};*/	/*64*/
/*FIELD2N poly_prime = {0x80000000,0x00000003};*/ /*63*/
/*FIELD2N poly_prime = {0x00000002,0x00002001};*/ /*33*/
/*FIELD2N poly_prime = {0x00000001,0x000000c5};*/ /*32*/
/*FIELD2N poly_prime = {0x80000009};*/	/*31*/
/*FIELD2N poly_prime = {0x800021};*/	/* 23*/
/*FIELD2N poly_prime = {0x100009};*/	/* 20*/
/*FIELD2N poly_prime = {0x20021};*/		/* 17*/
/*FIELD2N poly_prime = {0x1002B};*/		/* 16*/
/*FIELD2N poly_prime = {0x8003};*/		/* 15*/
/*FIELD2N poly_prime = {0x89};*/		/* 7*/
/*FIELD2N poly_prime = {0x43};*/		/*6*/
/*FIELD2N poly_prime = {0x29};*/		/*5*/
//FIELD2N poly_prime = {0x13};		/*4*/
//FIELD2N poly_prime = {0x7}; /*3*/
//------------------------------------------------------------------
void poly_rot_left(FIELD2N* a)
{
        INDEX i;
        ELEMENT bit,temp;

        bit = (a->e[0] & UPRBIT) ? 1L : 0L;
        for (i=NUMWORD; i>=0; i--) {
           temp = (a->e[i] & MSB) ? 1L : 0L;
           a->e[i] = ( a->e[i] << 1) | bit;
           bit = temp;
        }
        a->e[0] &= UPRMASK;
}
//------------------------------------------------------------------
void poly_rot_right(FIELD2N* a)
{
        INDEX i;
        ELEMENT bit,temp;

        bit = (a->e[NUMWORD] & 1) ? UPRBIT : 0L;
        SUMLOOP(i) {
           temp = ( a->e[i] >> 1)  | bit;
           bit = (a->e[i] & 1) ? MSB : 0L;
           a->e[i] = temp;
        }
        a->e[0] &= UPRMASK;
}
//------------------------------------------------------------------
void poly_mod_add(FIELD2N* a, FIELD2N* b, FIELD2N* c, FIELD2N* p)
{
	INDEX i;
	FIELD2N a_temp;
	FIELD2N b_temp;
	poly_copy(a, &a_temp);
	poly_copy(b, &b_temp);
	
	SUMLOOP(i) c->e[i] = a_temp.e[i] ^ b_temp.e[i];
}
//------------------------------------------------------------------
void poly_mod_sub(FIELD2N* a, FIELD2N* b, FIELD2N* c, FIELD2N* p)
{
	INDEX i;
	
	SUMLOOP(i) c->e[i] = a->e[i] ^ b->e[i];
}
//------------------------------------------------------------------
void poly_mul_shift(DBLFIELD* a)
{
	ELEMENT	*eptr, temp, bit;	/* eptr points to one ELEMENT at a time  */
	INDEX	i;
	
	eptr = &a->e[DBLWORD];		/* point to end, note: bigendian processing  */
	bit = 0;					/*  initial carry bit is clear */
	DBLLOOP (i)
	{
		temp = (*eptr << 1) | bit;		/* compute result as temporary  */
		bit = (*eptr & MSB) ? 1L : 0L;	/* get carry bit from shift  */
		*eptr-- = temp;					/* save new result  */
	}
}
//------------------------------------------------------------------
void poly_null(FIELD2N* a)
{
	INDEX i;
	
	SUMLOOP(i) a->e[i] = 0L;
}

void poly_dblnull(DBLFIELD* a)
{
	INDEX i;
	
	DBLLOOP(i) a->e[i] = 0L;
}
//------------------------------------------------------------------
void poly_copy(FIELD2N* from, FIELD2N* to)
{
	INDEX	i;
	
	SUMLOOP(i) to->e[i] = from->e[i];
}
//------------------------------------------------------------------
void sngltodbl(FIELD2N* from, DBLFIELD* to)
{
	INDEX i;
	
	poly_dblnull (to);
	SUMLOOP(i) to->e[DBLWORD - NUMWORD + i] = from->e[i];
}
//------------------------------------------------------------------
void dbltosngl(DBLFIELD* from, FIELD2N* to)
{
	INDEX i;
	
	SUMLOOP(i) to->e[i] = from->e[DBLWORD - NUMWORD + i];
}
//------------------------------------------------------------------
void poly_mul_partial(FIELD2N* a, FIELD2N* b, DBLFIELD* c)
{
	INDEX i, bit_count, word;
	ELEMENT mask;
	DBLFIELD B;
	
/*  clear all bits in result  */

	poly_dblnull(c);
	
/*  initialize local copy of b so we can shift it  */

	sngltodbl(b, &B);
	
/*  for every bit in 'a' that is set, add present B to result  */

	mask = 1;
	for (bit_count=0; bit_count<NUMBITS; bit_count++)
	{
		word = NUMWORD - bit_count/WORDSIZE;
		if (mask & a->e[word])
		{
			DBLLOOP(i) c->e[i] ^= B.e[i];
		}
		poly_mul_shift( &B);			/* multiply copy of b by x  */
		mask <<= 1;				/* shift mask bit up  */
		if (!mask) mask = 1;	/* when it goes to zero, reset to 1  */
	}
}
//------------------------------------------------------------------
void div_shift(DBLFIELD* a)
{
	ELEMENT	*eptr, temp, bit;
	INDEX	i;
	
	eptr = (ELEMENT*) &a->e[0];
	bit = 0;
	DBLLOOP (i)
	{
		temp = (*eptr>>1) | bit;		/*  same as shift left but  */
		bit = (*eptr & 1) ? MSB : 0L;	/*  carry bit goes off other end  */
		*eptr++ = temp;
	}
}
//------------------------------------------------------------------
INDEX log_2 (ELEMENT x)
{
	ELEMENT ebit, bitsave, bitmask;
	INDEX	k, lg2;
	
	lg2 = 0;
	bitsave = x;					/* grab bits we're interested in.  */
	k = WORDSIZE/2;					/* first see if msb is in top half */
	bitmask = -1L<<k;				/* of all bits  */
	while (k)
	{
		ebit = bitsave & bitmask;	/* did we hit a bit?  */
		if (ebit)					/* yes  */
		{
			lg2 += k;				/* increment degree by minimum possible offset */
			bitsave = ebit;			/* and zero out non useful bits  */
		}
	/*  The following two lines are slick and tricky.  We are binary searching,
		so k is divided by 2.  It's also the base amount we can add at any given
		search point.  The mask continuously searches for upper blocks being set.
		Once found, the remaining lower blocks are zeroed, to make sure we don't
		use bits which aren't the most significant.
	*/
		k /= 2;
		bitmask ^= (bitmask >> k);
	}
	return( lg2);
}
//------------------------------------------------------------------
INDEX degreeof(ELEMENT* t, INDEX dim)
{
	INDEX	degree, k;

/*  This is generic routine for arbitrary length arrays. use ELEMENT
	pointer to find first non-zero ELEMENT.  We will add degree from some base,
	so initial base is at least one WORDSIZE smaller than maximum possible.
*/	

	degree = dim * WORDSIZE;
	for (k=0; k<dim; k++)					/* search for a non-zero ELEMENT  */
	{
		if (*t) break;
		degree -= WORDSIZE;
		t++;
	}

/* for inversion routine, need to know when all bits zero.  */

	if (!*t) return(-1);		

/*  find most significant bit in this element */
	
	degree += log_2( *t);
	return(degree);
}
//------------------------------------------------------------------
void poly_div(DBLFIELD* top, FIELD2N* bottom, FIELD2N* quotient, FIELD2N* remainder)
{
	INDEX		deg_top;
	INDEX		deg_bot;
	INDEX		deg_quot;
	INDEX		bit_count;
	INDEX		i;
	INDEX		equot;
	ELEMENT		topbit;
	ELEMENT*	tptr;
	DBLFIELD	shift;
	
/*  Step 1: find degree of top and bottom polynomials.  */

	deg_top = degreeof( top, DBLWORD);
	deg_bot = degreeof( bottom, NUMWORD);
	
/*  prepare for null return and check for null quotient */

	poly_null (quotient);
	if (deg_top < deg_bot)
	{
		dbltosngl(top, remainder);
		return;
	}

/*  Step 2: shift bottom to align with top.  Note that there 
		are much more efficient ways to do this. */
	
	deg_quot = deg_top - deg_bot;
	bit_count = deg_quot + 1;
	sngltodbl(bottom, &shift);
	for (i = 0; i<deg_quot; i++)
		 poly_mul_shift(&shift);

/* Step 3: create bit mask to check if msb of top is set  */

	topbit = 1L << (deg_top % WORDSIZE);
	tptr = (ELEMENT*) top + DBLWORD - deg_top/WORDSIZE;
	
/*  for each possible quotient bit, see if bottom can be subtracted
	(added, it's modulo 2!) from top. If it can, set that bit in
	quotient. If it can't, clear that bit in quotient.  Shift bottom
	right once and keep going for total bit_count.
*/

	while (bit_count)
	{
	
/*  Step 4: determine one bit of the quotient.  */

		if (*tptr & topbit)		/* is bit set in top?  */
		{
			DBLLOOP (i)			/* yes, subtract shift from top  */
				top->e[i] ^= shift.e[i];
				
	/* find word and bit in quotient to be set  */
	
			equot = 	NUMWORD - deg_quot/WORDSIZE;	
			quotient->e[equot] |= 1L << (deg_quot % WORDSIZE);	
		}

/*  Step 5: advance to the next quotient bit and perform the necessary
			shifts.  */
			
		bit_count--;			/*  number of bits is one more  */
		deg_quot--;				/*  than polynomial degree  */
		div_shift(&shift);		/*  divide by 2  */
		topbit >>= 1;			/*  move mask over one bit also  */
		if (!topbit)
		{
			topbit = MSB;		/*  reset mask bit to next word  */
			tptr++;				/*  when it goes to zero  */
		}
	}

/*  Step 6: return the remainder in FIELD2N size  */

	dbltosngl (top, remainder);
}
//------------------------------------------------------------------
void poly_mod_mul(FIELD2N *a, FIELD2N* b, FIELD2N* c, FIELD2N* p)
{
	DBLFIELD temp;
	FIELD2N  dummy;
	
	poly_mul_partial(a, b, &temp);
	poly_div(&temp, p, &dummy, c);
}
//------------------------------------------------------------------
void poly_mod_inv(FIELD2N* a, FIELD2N* inverse, FIELD2N* p)
{
	FIELD2N	pk, pk1, pk2;
	FIELD2N rk, rk1;
	FIELD2N qk, qk1, qk2;
	INDEX	i;
	DBLFIELD rk2;
	
/*  initialize remainder, quotient and product terms  */

	sngltodbl (p, &rk2);	
	poly_copy( a, &rk1);
	poly_null( &pk2);
	poly_null( &pk1);
	pk1.e[NUMWORD] = 1L;
	poly_null( &qk2);
	poly_null( &qk1);
	qk1.e[NUMWORD] = 1L;

/*  compute quotient and remainder for Euclid's algorithm.
    when degree of remainder is < 0, there is no remainder, and we're done.
    At that point, pk is the answer.
*/
	poly_null( &pk);
	pk.e[NUMWORD] = 1L;
	poly_div(&rk2, &rk1, &qk, &rk);
		
	while (degreeof(&rk, NUMWORD) >= 0)
	{
		poly_mul_partial( &qk, &pk1, &rk2);
		SUMLOOP(i) pk.e[i] = rk2.e[i+DBLWORD-NUMWORD] ^ pk2.e[i];
		
/*  set up variables for next loop */

		sngltodbl(&rk1, &rk2);
		poly_copy( &rk, &rk1);
		poly_copy( &qk1, &qk2);
		poly_copy( &qk, &qk1);
		poly_copy( &pk1, &pk2);
		poly_copy( &pk, &pk1);
		poly_div( &rk2, &rk1, &qk, &rk);
	}
	poly_copy( &pk, inverse);		/* copy answer to output  */
}
//------------------------------------------------------------------
void poly_gcd(FIELD2N* u, FIELD2N* v, FIELD2N* gcd)
{
	DBLFIELD	top;
	FIELD2N		r, dummy, temp;
	
	sngltodbl( u, &top);
	poly_copy( v, &r);
	
	while( degreeof( &r, NUMWORD) >= 0)
	{
		poly_div( &top, &r, &dummy, &temp);
		sngltodbl( &r, &top);
		poly_copy( &temp, &r);
	}
	dbltosngl( &top, gcd);
}
//------------------------------------------------------------------
void mul_x_mod(FIELD2N* u, FIELD2N* v, FIELD2N* w)
{
	DBLFIELD mulx;
	INDEX    i, deg_v;
	
	deg_v = degreeof(v, NUMWORD);
	sngltodbl( u, &mulx);
	poly_mul_shift( &mulx);  /*  multiply u by x  */
	dbltosngl( &mulx, w);
	if (w->e[ NUMWORD - (deg_v/WORDSIZE)] & ( 1L << (deg_v % WORDSIZE)))
		SUMLOOP (i) w->e[i] ^= v->e[i];
}
//------------------------------------------------------------------
INDEX irreducible(FIELD2N* v)
{
	FIELD2N		vprm, gcd, x2r, x2rx, temp;
	FIELD2N		sqr_x[NUMBITS+1];
	INDEX		i, r, deg_v, k;
		
/*  check that gcd(v, v') = 1.  If not, then v not irreducible.  */ 

	SUMLOOP(i) vprm.e[i] = (v->e[i] >> 1) & DERIVMASK;
	poly_gcd( v, &vprm, &gcd);
	if (gcd.e[NUMWORD] > 1) return(0);
	for (i=0; i<NUMWORD; i++)  if (gcd.e[i]) return(0);

/*  find maximum power we have to deal with  */

	deg_v = degreeof(v, NUMWORD);

/*  create a vector table of powers of x^2k mod v.
	this will be used to compute square of x^2^r
*/

	poly_null (&sqr_x[0]);
	sqr_x[0].e[NUMWORD] = 1;
	for (i=1; i<= deg_v; i++)
	{
		mul_x_mod( &sqr_x[i-1], v, &temp);
		mul_x_mod( &temp, v, &sqr_x[i]);
	}

/*  check that gcd( x^2^r - x, v) == 1 for all r <= degreeof(v)/2.
	set x^2^r = x for r = 0 to initialize.
*/
	poly_null( &x2r);
	x2r.e[NUMWORD] = 2;
	for ( r=1; r <= deg_v/2; r++)
	{
/*  square x^2^r mod v.  We do this by seeing that s^2(x) = s(x^2).
	for each bit j set in x^2^(r-1) add in x^2j mod v to find x^2^r.
*/
		poly_null( &x2rx);
		for (i=0; i <= deg_v; i++)
		{
			if ( x2r.e[NUMWORD - (i/WORDSIZE)] & ( 1L << (i%WORDSIZE) ))
				SUMLOOP(k) x2rx.e[k] ^= sqr_x[i].e[k];
		}
		
/*  save new value of x^2^r mod v and compute x^2^r - x mod v */

		poly_copy( &x2rx, &x2r);
		x2rx.e[NUMWORD] ^= 2;

/*  is gcd( x^2^r - x, v) == 1?  if not, exit with negative result  */

		poly_gcd( &x2rx, v, &gcd);
		if (gcd.e[NUMWORD] > 1) return (0);
		for (i=0; i<NUMWORD; i++) if ( gcd.e[i]) return (0);
	}

/*  passed all tests, provably irreducible  */

	return (1);
}
//------------------------------------------------------------------
void poly_mon_mod_mul(FIELD2N* a, FIELD2N* b, FIELD2N* c, FIELD2N* p)
{
	INDEX   i;
	FIELD2N two;
	FIELD2N two_inv;
	FIELD2N a_temp;
	FIELD2N b_temp;

	poly_copy(a, &a_temp);
	poly_copy(b, &b_temp);

	poly_null(&two);
	two.e[NUMWORD] = 0x00000002;

	poly_mod_inv(&two, &two_inv, p);

	poly_mod_mul(&a_temp, &two_inv, c, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		poly_mod_mul(c, &two_inv, c, p);
	}
	poly_mod_mul(c, &b_temp, c, p);
}
//------------------------------------------------------------------
void poly_mon_mod_div(FIELD2N* a, FIELD2N* b, FIELD2N* c, FIELD2N* p)
{
	INDEX   i;
	FIELD2N two;
	FIELD2N b_inv;
	FIELD2N temp;

	poly_null(&two);
	two.e[NUMWORD] = 0x00000002;

	poly_mod_inv(b, &b_inv, p);
	poly_mod_mul(a, &b_inv, c, p);
	poly_mod_mul(c, &two, c, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		poly_mod_mul(c, &two, &temp, p);
		poly_copy(&temp, c);
	}
}
//------------------------------------------------------------------
void poly_to_mon(FIELD2N* a, FIELD2N* b, FIELD2N* p)
{
	INDEX   i;
	FIELD2N two;

	poly_null(&two);
	two.e[NUMWORD] = 0x00000002;

	poly_mod_mul(a, &two, b, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		poly_mod_mul(b, &two, b, p);
	}
}
//------------------------------------------------------------------
void mon_to_poly(FIELD2N* a, FIELD2N* b, FIELD2N* p)
{
	INDEX   i;
	FIELD2N two;
	FIELD2N two_inv;

	poly_null(&two);
	two.e[NUMWORD] = 0x00000002;

	poly_mod_inv(&two, &two_inv, p);

	poly_mod_mul(a, &two_inv, b, p);
	for(i = 0; i < NUMBITS-1; i++)
	{
		poly_mod_mul(b, &two_inv, b, p);
	}
}
//------------------------------------------------------------------
void print_poly(FIELD2N* a)
{
	INDEX i;
	
	SUMLOOP(i) printf("%8x ", a->e[i]);
	printf("\n");
}
//------------------------------------------------------------------
void wfile_poly(FIELD2N* a, char* filename)
{
	FILE* cfPtr;
	INDEX i;

	if((cfPtr = fopen(filename, "w")) == NULL) //"a" => append
	{
		printf("File could not be opened\n");
	}
	else
	{
		printf("Write %s file\n", filename);
		SUMLOOP(i) fprintf(cfPtr, "%08x", a->e[i]);
	}
	fprintf(cfPtr, "\n");

	fclose(cfPtr);
}
//------------------------------------------------------------------
void afile_poly(FIELD2N* a, char* filename)
{
	FILE* cfPtr;
	INDEX i;

	if((cfPtr = fopen(filename, "a")) == NULL) //"a" => append
	{
		printf("File could not be opened\n");
	}
	else
	{
		printf("Write %s file\n", filename);
		SUMLOOP(i) fprintf(cfPtr, "%08x", a->e[i]);
	}
	fprintf(cfPtr, "\n");

	fclose(cfPtr);
}
